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We present an atomic orbital based approximate scheme for self-interaction correction (SIC) to 
the local density approximation of density functional theory. The method, based on the idea of 
Filippetti and Spaldin [Q Phys. Rev. B 67, 125109 (2003)], is implemented in a code using 
localized numerical atomic orbital basis sets and is now suitable for both molecules and extended 
solids. After deriving the fundamental equations as a non-variational approximation of the self- 
consistent SIC theory, we present results for a wide range of molecules and insulators. In particular, 
we investigate the effect of re-scaling the self-interaction correction and we establish a link with 
the existing atomic-like corrective scheme LDA+[/. We find that when no re-scaling is applied, 
i.e. when we consider the entire atomic correction, the Kohn-Sham HOMO eigenvalue is a rather 
good approximation to the experimental ionization potential for molecules. Similarly the HOMO 
eigenvalues of negatively charged molecules reproduce closely the molecular affinities. In contrast a 
re-scaling of about 50% is necessary to reproduce insulator bandgaps in solids, which otherwise are 
largely overestimated. The method therefore represents a Kohn-Sham based single-particle theory 
and offers good prospects for applications where the actual position of the Kohn-Sham eigenvalues 
is important, such as quantum transport. 

PACS numbers: 



I. INTRODUCTION 

Density functional theory (DFT), in both its static Q 
and time-dependent Q forms, has become by far the 
most successful and widely used among all the electronic 
structure methods. The most obvious reason for this suc- 
cess is that it provides accurate predictions of numerous 
properties of atoms, inorganic molecules, bio-molecules, 
nanostructures and solids, thus serving different scientific 
communities. 

In addition DFT has a solid theoretical foundation. 
The Hohenberg-Kohn theorem |2( establishes the exis- 
tence of a unique energy functional E[p] of the electron 
density p which alone is sufficient to determine the ex- 
act ground-state of a TV-electron system. Although the 
energy functional itself is not known, several of its gen- 
eral properties can be demonstrated rigorously. These 
are crucial for constructing increasingly more predictive 
approximations to the functional and for addressing the 
failures of such approximations Q. 

Finally, but no less important, the Kohn-Sham (KS) 
formulation of DFT |f| establishes a one to one map- 
ping of the intrinsically many-body problem onto a ficti- 
tious single-particle system and offers a convenient way 
for minimizing E[p\. The degree of complexity of the 
Kohn-Sham (KS) problem depends on the approxima- 
tion chosen for the density functional. In the case of the 
local density approximation (LDA) [j| the KS problem 
typically scales as TV 3 , where the scaling is dominated by 
the diagonalization algorithm. However, clever choices 
with regards to basis sets and sophisticated numerical 
methods make order- TV scaling a reality 0, Q • 

The energy functional E[p^ , p^] (p a , a =T, J. is the spin 



density, p = J2 a P a ) can be written as 
E[p\p±}=T s [p] + J d 3 rp(r)v(r) + U{p} + E xc [p\pt}, 

w 

where Tg is the kinetic energy of the non-interacting sys- 
tem, v(r) the external potential, U the Hartree electro- 
static energy and E xc the exchange and correlation (XC) 
energy. This last term is unknown and must be approxi- 
mated. The construction of an approximated functional 
follows two strategies: empirical and "constraint satis- 
faction" . 

Empirical XC functionals usually violate some of the 
constraints imposed by exact DFT, and rely on parame- 
tcrizations obtained by fitting representative data. One 
includes in this category, XC functionals which borrow 
some functional dependence from other theories. This is 
for instance the case of the celebrated LDA+U scheme 
0, 0, where the Hubbard-t/ energy takes the place of 
the LDA energy for certain "strongly correlated" atomic 
orbitals (typically d and / shells). The method however 
depends on the knowledge of the Coulomb and exchange 
parameters U and J, which vary from material to mate- 
rial, and can also be different for the same ion in different 
chemical environments 0, . 

In contrast the construction based on "constraint sat- 
isfaction" proceeds by developing increasingly more so- 
phisticated functionals, which nevertheless satisfy most 
of the properties of exact DFT ^1 • It was argued that 
this method constructs a "Jacob's ladder" [LJ, where 
functionals are assigned to different rungs depending 
on the number of ingredients they include. Thus the 
LDA, which depends only on the spin-densities is on 
the first rung, the generalized gradient approximation 
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(GGA) 14|, which depends also on Vp ff , is on the second 
rung, the so-called meta-GGA functionals which in 
addition to p a and V ' p® depend on either the Laplacian 
'V 2 p a or the orbital kinetic energy density, are on the 
third rung and so on. The higher its position on the 
ladder the more accurate a functional becomes, but at 
the price of increasing computational overheads. There- 
fore its worth investigating corrections to the functionals 
of the lower rungs, which preserve most of the funda- 
mental properties of DFT and do not generate massive 
additional numerical overheads. 

One of the fundamental problems intrinsic to the semi- 
local functionals of the first three rungs is the presence of 
self- interaction (SI) [Ifij . This is the spurious interaction 
of an electron in a given KS orbital with the Hartree and 
XC potential generated by itself. Such an interaction 
is not present in the Hartree-Fock method, where the 
Coulomb self-interaction of occupied orbitals is cancelled 
exactly by the non-local exchange. However when using 
semi-local functionals such a cancellation is not complete 
and the rigorous condition for KS-DFT 



U{p°}+ E xc [p°, 0} = Q 



(2) 



for the orbital density p^ = \ipn\ 2 °f the fully occupied 
KS orbital is not satisfied. A direct consequence of 
the self- interaction in LDA/GGA is that the KS potential 
becomes too repulsive and exhibits an incorrect asymp- 
totic behavior [lq . 

This "schizophrenic" (self-interacting) nature of semi- 
local KS potentials generates a number of failures in de- 
scribing elementary properties of atoms, molecules and 
solids. Negatively charged ions (H~, 0~, F~) and 
molecules are predicted to be unstable by LDA [ljfl and 
transition metal oxides are described as small-gap Mott- 
Hubbard antiferromagnets (MnO, NiO) [3 

or even as 

ferromagnetic metals (FeO, CoO) [18J instead of charge- 
transfer insulators. Moreover the KS highest occupied 
molecular orbital (HOMO), the only KS eigenvalue that 
can be rigorously associated to a single particle energy 
[l^. I20L I21J. is usually nowhere near the actual ionization 
potential |l6j. 

Finally XC functionals affected by SI do not present 
a derivative discontinuity as a function of the occupa- 
tion [T^.l20| . Semi- local functionals in fact continuously 
connect the orbital levels of systems of different integer 
occupation. This means for instance that when adding 
an extra electron to an open shell TV-electron system the 
KS potential does not jump discontinuously by In — 
where In and An are respectively the ionization potential 
and the electron affinity for the ./V-electron system. This 
serious drawback is responsible for the incorrect disso- 
ciation of heteronuclear molecules into charged ions [2^ 
and for the metallic conductance of insulating molecules 

The problem of removing the SI from a semi-local po- 
tential was acknowledged a long time ago when Fermi and 
Amaldi proposed a first crude correction [24[. However 
the modern theory of self-interaction corrections (SIC) in 



DFT is due to the original work of Perdew and Zunger 
from almost three decades ago 0]. Their idea consists 
in removing directly the self-Hartree and self-XC energy 
of all the occupied KS orbitals from the LDA XC func- 
tional (the same argument is valid for other semi-local 
functionals), thus defining the SIC functional as 

occupied 

(3) 

Although apparently simple, the SIC method is more 
involved than standard KS DFT. The theory is still a 
density functional one, i.e. it satisfies the Hohenberg- 
Kohn theorem, however it does not fit into the Kohn- 
Sham scheme, since the one-particle potential is orbital- 
dependent. This means that one cannot define a kinetic 
energy functional independently from the choice of E xc 
|l6j| . Two immediate consequences are that the ip° are 
not orthogonal, and that the orbital-dependent potential 
can break the symmetry of the system. This last aspect 
is particularly important for solids since one has to give 
up the Bloch representation. 

In this paper we explore an approximate method for 
SIC to the LDA, which has the benefit of preserving the 
local nature of the LDA potential, and therefore main- 
tains all of the system's symmetries. We have followed in 
the footsteps of Filippetti and Spaldin pj , who extended 
the original idea of Vogel and co-workers [2M l26l 
of considering only the atomic contributions to the SIC. 
We have implemented such a scheme into the localized 
atomic orbital code Siesta [2^] and applied it to a vast 
range of molecules and solids. In particular we have in- 
vestigated in detail how the scheme performs as a single- 
particle theory, and how the SIC should be rescaled in 
different chemical environments. 



II. REVIEW OF EXISTING METHODS 

The direct subtraction proposed by Perdew and Zunger 
is the foundation of the modern SIC method. However 
the minimization of the SIC functional is not trivial, 
in particular for extended systems. The main problem 
is that E xc itself depends on the KS orbitals. Thus it 
does not fit into the standard KS scheme and a more 
complicated minimization procedure is needed. 

Following the minimization strategy proposed by Levy 
|29|. which prescribes to minimize the functional first 
with respect to the KS orbitals and then with re- 
spect to the occupation numbers p^, Perdew and Zunger 
derived a set of single-particle equations 
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<ff,nO) 



Ml = <r,SIC , a 



(4) 



where the effective single-particle potential v° s n (r) is de- 
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fined as 



,ct,LDA 



([p\p l ] 



,(r) = v(r) +u([p];r) - 

- U ([p„];r)-<f DA ([pT )0 ];r), (5) 



and 



*(W;r) = | d 3 r' 



|r-r' 



<5,o ff (r)' 



(0) 
(7) 



These are solved in the standard KS way for atoms, 
with good results in terms of quasi-particle energies ^(| . 
In this particular case the KS orbitals ip° show only small 
deviations from orthogonality, which is re-imposed with 
a standard Schmidt orthogonalization. 

The problem of the non-orthogonality of the KS or- 
bitals can be easily solved by imposing the orthogonality 
condition as a constraint to the energy functional, thus 
obtaining the following single-particle equation 



1 

-2 V ' 



t (r) 



(8) 



The first consists in a direct minimization with respect 
to the localized orbitals <f%, i.e. in solving equation © 
when we replace tp with (f> and the orbital densities en- 
tering the definition of the one-particle potential JSJ are 
simply Pn — \4>n\ 2 ■ In addition the following minimiza- 
tion condition must be satisfied 



„SiC 



v 



SIC IC 



) = o 



(9) 



where = u([p n ];r) + v^ DA ([p} l , 0]; r). An expres- 
sion for the gradient of the SIC functional, which also 
constrains the orbitals to be orthogonal has been derived 
[33| and applied to atoms and molecules with a mixture 
of successes and bad failures [3 HE EH • 

The second strategy uses the canonical orbitals i> and 
seeks the minimization of the SIC energy by varying both 
the orbitals if) and the unitary transformation M.. The 
corresponding set of equations is 



M nm (j>° 



(10) 



(11) 



Even in this case where orthogonality is imposed, two 
major problems remain: the orbitals minimizing the en- 
ergy functional are not KS-type and in general do not 
satisfy the system's symmetries. 

If one insists in minimizing the energy functional in a 
KS fashion by constructing the orbitals according to the 
symmetries of the system, the theory will become size- 
inconsistent, or in other words it will be dependent on 
the particular representation employed. Thus one might 
arrive at a paradox, where in the self-interaction of N 
hydrogen atoms arranged on a regular lattice of large 
lattice spacing (in such a way that there is no interaction 
between the atoms) vanishes, since the SIC of a Bloch 
state vanishes for N — > oo. However the SIC for an 
individual H atom, when calculated using atomic-like or- 
bitalSj_accounts for essentially all the ground-state energy 
error [lfj . Therefore a size-consistent theory of SIC DFT 
must look for a scheme where a unitary transformation of 
the occupied orbitals, which minimizes the SIC energy is 
performed. This idea is at the foundation of all modern 
implementations of SIC. 

Significant progress towards the construction of a size- 
consistent SIC theory was made by Pederson, Heaton 
and Lin, who introduced two sets of orbitals: local- 
ized orbitals <^ minimizing E^P and canonical (Kohn- 
Sham) de-localized orbitals [H HU, HI . The local- 
ized orbitals are used for defining the densities entering 
into the effective potential ©, while the canonical or- 
bitals are used for extracting the Lagrangian multipli- 
ers e^ IC , which are then associated to the KS eigenval- 
ues. The two sets are related by unitary transformation 
= Em-^nmCi an( l one nas two possible strategies 
for minimizing the total energy. 



SIC fm 



(12) 



where Hq is the standard LDA Hamiltonian (without 
SIC). Thus the SIC potential for the canonical orbitals 
appears as a weighted average of the SIC potential for the 
localized orbitals. The solutions of the set of equations 
(|10|l is somehow more appealing than that associated to 
the localized orbitals since the canonical orbitals can be 
constructed in a way that preserves the system's symme- 
tries (for instance translational invariance). 

A convenient way for solving the equation IjlQI) is that 
of using the so called "unified Hamiltonian" method [3(| . 
This is defined as (we drop the spin index a) 



occup 



P n H^P n - 



occup 

E 



(P n H n Q + QH n P n ) + QHqQ , 

(13) 

where P n = IVO^nl ^ s the projector over the occu- 
pied orbital ip^i an( l Q i s l ne projector over the unoc- 
cupied ones Q — 1 — X)°i° CUP Pn- The crucial point is that 
the diagonal elements of the matrix e^ IC and their cor- 
responding orbitals are respectively eigenvalues and 
eigenvectors of H u , from which the whole e^ IC can l> e 
constructed. Finally, and perhaps most important, at the 
minimum of the SIC functional, the canonical orbitals di- 
agonalize the matrix e^ IC , whose eigenvalues can now 
be interpreted as an analogue of the Kohn-Sham eigen- 
values [32| . 

It is also interesting to note that an alternative way for 
obtaining orbital energies is that of constructing an ef- 
fective Sl-free local potential using the Krieger-Li-Iafrate 
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method l37l . This has been recently explored by several 
groups |38l l39l Eof 

When applied to extended systems the SIC method de- 
mands considerable additional computational overheads 
over standard LDA. Thus for a long time it has not en- 
countered the favor of the general solid state community. 
In the case of solids the price to pay for not using canon- 
ical orbitals is enormous since the Bloch representation 
should be abandoned and in principle infinite cells should 
be considered. For this reason the second minimization 
scheme, in which the canonical orbitals are in a Bloch 
form, is more suitable. In this case for each k-vector 
one can derive an equation identical to equation (|l(Jfl. 
where e^ IC = e^ IC (k) and n is simply the band index 
[4l|. The associated localized orbitals cf> for instance can 
be constructed as Wannier functions and the minimiza- 
tion scheme proceeds in a similar way to that done for 
molecules. 

The problem here is that in practice, the cell needed 
to describe the localized states </> may be considerably 
larger than the primitive unit cell. This is not the case 
for ionic insulators where the localized orbitals are 
well approximated by atomic orbitals. Such a simplifi- 
cation is however not valid in general. For example su- 
percells as large as 500 atoms have been considered for 
describing the localized d shells in transition metals ox- 
ides 0>EJS|- Despite these difficulties the SIC scheme 
has been applied to a vast range of solid state systems 
with systematic improvement over LDA. These include, 
in addition to transition metals monoxides |42l 0, , 
rare-earth materials diluted magnetic semiconduc- 
tors F e 304 0, heavy elements compounds |49j. iust 
to name a few. 

In order to make the SIC method more scalable sev- 
eral approximations have been proposed. One possibility 
is that of incorporating part of the SIC into the defi- 
nition of the pseudopotentials The idea consists 
in subtracting the atomic SI from the free atom, and 
then transferring the resulting electronic structure to the 
definition of a standard norm-conserving pseudopoten- 
tial. This approximation is sustained by the fact that 
the transformation 7W, which relates canonical and lo- 
calized orbitals does hardly mix core and valence states 
[32|. Thus the SIC contribution to the total energy can 
be separated into the contributions from the core and 
the valence and in first approximation the latter can be 
neglected |5l|. The benefit of this method is that trans- 
lational invariance is regained and the complicated pro- 
cedure of minimizing M. is replaced by a pseudopotential 
calculation. 

A further improvement over the pseudopotential 
method was recently presented by Vogel and co-workers 
[25L 12a . [23 and then extended by Filippetti and Spaldin 
Q. The method still assumes separability between the 
core and the valence contributions to the SIC, but now 
the SIC for the valence electrons is approximated by an 
atomic- like contribution, instead of being neglected. This 
atomic SIC (ASIC) scheme is certainly a drastic approx- 



imation, since it implicitly assumes that the transforma- 
tion M. minimizing the SIC functional leads to atomic 
like orbitals. 

In the work of Vogel this additional contribution is not 
evaluated self-consistently for the solid, while the imple- 
mentation of Filippetti assumes a linear dependance of 
the SIC over the orbital occupation. In spite of the ap- 
proximations involved, the method has been applied suc- 
cessfully to a range of solids including II- VI semiconduc- 
tors and nitrites [J, l25l |26l , transition metal monoxides 
P, , silver halides T27l . noble metal oxides fer- 
roelectric materials [H l54l 155) . high-k materials |56l| and 
diluted magnetic semiconductors • Interestingly 

most of the systems addressed by the ASIC method are 
characterized by semi-core d orbitals, for which an atomic 
correction looks appropriate, and a similar argument is 
probably valid for ionic compounds as recently demon- 
strated for the case of SiC £3 ■ 

Here we further investigate the self-consistent ASIC 
method PJ by examining both finite and extended sys- 
tems, and by critically considering whether a scaling fac- 
tor, additional to the orbital occupation, is needed for 
reproducing the correct single particle spectrum. 



III. FORMALISM AND IMPLEMENTATION 

In this section we derive the fundamental equations 
of the ASIC method, while looking closely at the main 
approximations involved in comparison to the fully self- 
consistent SIC approach. Our practical implementation 
is also described. 



A. The ASIC potential 

The starting point of our analysis is the SIC 
Schrodinger-like equation (fTUf) for the canonical orbitals. 
Let us assume, as from reference [HlJ, that the rotation 
M. transforming localized orbitals (to be determined) 
into canonical orbitals (see equation Ulljl) does not mix 
core and valence states. We also assume that core elec- 
trons are well localized into atomic-like wave-functions 
and that they can be effectively described by a norm- 
conserving pseudopotential. 

Let us now assume that M. is known and so are the 
localized orbitals ^m- In this case the canonical orbitals 
diagonalize e^ IC and the equation ifTUl) simply reduces 
to 

(ftf + At£ IC )< = # SIC <, (14) 

with Au„ defined in equation i|12|l . The Hamiltonian 
Hq + Aw^ IC can be then re-written in a convenient form 

as 

Hg + Av^^HS + ^v^P*, (15) 

m 
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where vf^ is the self-interaction potential for the local- 
ized orbital and P^ is the projector over the same 
state 

i*<(r) = CM / dV<(r')CV) = C(r)<CIO ■ 

(16) 

Two main approximations are then taken in the ASIC 
approach 0, [23. First the localized states § a m are as- 
sumed to be atomic-like orbitals (ASIC orbitals) and 
the SIC potential is approximated as 

£t£ c (r)i» ^E^. ( 17 ) 

m m 

with < src (r) = u([p m ];r) + v^ DA ([pl,0];v) and = 
|$^| 2 , P* is the projector of equation ((TBj) obtained by 
replacing the </>'s with the ASIC orbitals and a is a 
scaling factor. Note that the orbitals & m are not explic- 
itly spin-dependent and one simply has $^=$ m p^j with 
■p a m the orbital occupation (p^ = 0, 1). The factor a is an 
empirical factor, which accounts for the particular choice 
of ASIC orbitals. This first approximation is expected to 
be accurate for systems retaining an atomic-like charge 
density as in the case of small molecules. It is also for- 
mally exact in the one-electron limit (for a = 1). In 
the case of extended solids the situation is less transpar- 
ent, since in general the functions minimizing are 
Wannier-like functions 60] . 

The second approximation taken in the ASIC method 
is that of replacing the non-local projector P* with its 
expectation value. The idea is that the SIC potential for 
the canonical orbitals Aw^ IC is formally a weighted av- 
erage of the SIC potential for the localized orbitals wf c . 
For the exact SIC method the weighting factor is the 
non-local projector Mnm^fir- This means that the SIC 
potential for a given canonical orbital iji n is maximized in 
those regions where the overlap with some of the localized 
orbitals <p m is maximum. In the ASIC method such non- 
local projector is replaced more conveniently by a scalar. 
In the original proposal by Voegl et al. [251 12a W\ this 
was simply set to one. Here we consider the orbital occu- 
pation pf n of the given ASIC orbital $ m , i.e. we replace 
P* with its expectation value 

pi - =f m =£#«i^io , (18) 

n 

where f% is the occupation number of the Kohn-Sham 
orbital 1/;%,. The final form of ASIC potential is then 

«Wr) = «£< SIC ( r K>' (19) 

m 

Let us now comment on the empirical scaling factor 
a. In reference 0] a was set to 1/2 in order to capture 
eigenvalue relaxation. This choice however violates the 
one-electron limit of the SIC potential, which is correctly 
reproduced for a = 1. We can then interpret a as a 



measure of the deviation of the ASIC potential from the 
exact SIC potential. Ultimately a reflects the deviation 
of the actual ASIC projectors |3>)(3>| from the localized 
orbitals defining the SI corrected ground state. One then 
expects a to be close to 1 for systems with an atomic-like 
charge density, and to vanish for metals, whose valence 
charge density resembles that of a uniform electron gas 
|6lj | . Intermediate values are then expected for situations 
different from these two extremes, and we will show that 
a values around 1/2 describe well a vast class of mid- to 
wide-gap insulators. 



B. Implementation 

The final form of the SIC potential to subtract from 
the LSDA (local spin density approximation) one (equa- 
tion 119|) ) is that of a linear combination of non-local 
pseudopotential-like terms. These are uniquely defined 
by the choice of exchange and correlation potential used 
and by ASIC orbitals $ m . The practical way of con- 
structing such potentials, i.e. the way of importing the 
atomic SIC to the solid state, depends on the specific nu- 
merical implementation used for the DFT algorithm. At 
present only plane-wave implementations are available 
IH I25L Hi| [23, while here we present our new scheme 
based on the pseudo atomic orbital (PAO) 65] code 
Siesta $2f£. 

We start by solving the atomic all-electron SIC-LSDA 
equation for all the species involved in the solid state 
calculation. Here we apply the original Perdew-Zunger 
(PZ-SIC) formalism [l(j and we neglect the small non- 
orthogonality between the Kohn-Sham orbitals. Thus 
we obtain a set of SI corrected atomic orbitals $ m , which 
exactly solve the atomic SIC-LSDA problem. The atomic 
orbitals <& m describing the valence electrons are then used 
to define the ASIC potentials v^ IC 

<f IC (r) = u([p m ];r) + v^ UA ([pl, 0]; r) (20) 

with P ° m = |$ m | 2 . 

At the same time a standard LSDA calculation for the 
same atoms is used to construct the pseudopotentials de- 
scribing the core electrons. These are standard norm- 
conserving scalar relativistic Troullier-Martins pseudopo- 
tentials |62j ] with nonlinear core corrections |63| . Thus we 
usually neglect the SIC over the core states, when con- 
structing the pseudopotentials. This is justified by the 
fact that the eigenvalues for the SIC-LSDA-pseudoatom, 
i.e. for the free atom where the effects of core electrons 
are replaced by LSDA pseudopotentials but SIC is ap- 
plied to the valence electrons, are in excellent agreement 
with those obtained by all-electron SIC-LSDA calcula- 
tions m. 

The final step is that of recasting the ASIC potentials 
"D^j C (r), which have a — 2/r asymptotic behaviour, in a 
suitable non-local form. This is obtained with the stan- 
dard Kleinman-Bylander [64| scheme and the final ASIC 



6 



potential (equation JT5Jl) is written as 

\l°m)(l a m \ 



V ASIC — X] ' 



no 



where the ASIC projectors are given by 

7 -(r)=a^C SIC (r)^(r) 
and the normalization factors are 



Co 



ap ri 



i^sic 



(21) 



(22) 



(23) 



The orbital functions $ m are atomic-like functions 
with a finite range, which ensure that the ASIC projec- 
tors 7 m do not extend beyond that range. These are con- 
structed in the same way as the Siesta basis set orbitals, 
i.e. as solutions of the pseudo-atomic problem with an 
additional confining potential at the cutoff radius r cuto ff 
[65|. The choice of the appropriate cutoff radius for the 
SIC projectors should take into account the two following 
requirements. On the one hand it should be sufficiently 
large to capture most of the SIC corrections. A good 
criterion p| is that the SIC-LSDA contribution to the 
orbital energy of the free atom 



SIC m 



7Tl/ 



(24) 



is reproduced within some tolerance. On the other hand 
the cutoff should be reasonably short so as not to change 
the connectivity of the matrix elements of the PAO 
Hamiltonian. In other words we need to ensure that or- 
bitals otherwise considered as disconnected in evaluating 
the various parts of the Hamiltonian matrix are not con- 
sidered connected for the «asic ma trix elements alone. 

As a practical rule we set the cutoff radius for a par- 
ticular orbital of a given atom to be either equal to the 
longest among the cutoff radii of the PAO basis set for 
that particular atom (typically the first £ of the lowest 
angular momentum), or, if shorter, the radius at which 
fegic m < O.lmRy. Typically, when reasonable cutoff 
radii (6 to 9 Bohr) are used, we find that the atomic 
SIC-LSDA eigenvalues are reproduced to within 1 to 5 
mRy for the most extended shells and to within 0.1 mRy 
for more confined shells. Thus 5s'. 



SIC m 



are rather well 
converged already for cutoff radii defined by a PAO ener- 
gies shifts [28| of around 20mRy, although usually smaller 
PAO energy shifts are necessary for highly converged to- 
tal energy calculations. 

Finally the matrix elements of the SIC potential are 
calculated as usual over the Siesta basis set. Additional 
basis functions Xm are constructed from the confined lo- 
calized atomic orbitals described before using the split- 
norm scheme |65l I66L l67| . The density matrix p a is rep- 
resented over such basis p" and the orbital populations 
are calculated as 



Pm ^ ] S m fj,Pp U S v 



(25) 



where S mfi is a matrix element of the overlap matrix. 
Note that in principle the orbital population should be 
constructed for the atomic SIC orbital <E> m . However, we 
notice that p m is rather insensitive to the specific choice 
of orbital, once this has a reasonable radial range. For 
practical numerical reasons in the present implementa- 
tion, we always use the orbital populations projected onto 
the basis set sub-space consisting of the most extended 
first-£ orbitals of the atomic species involved. The matrix 
elements of the SIC potential are simply 



( w ASIc) 



(X/x|7m)(7mlX*) 

Cm 



(26) 



and the ASIC-KS equation takes the final form 



1 



,LSDA 



V ASIC 



lb" = t< J ' SlC ib' T 



with upp the pseudopotcntial. 



(27) 



C. Total Energy 

The ener gy c orresponding to the SIC-LSDA functional 
is given by [lj| 

occ. 

(28) 



where 



U[Pn] = / d 3 r -p n {v)u{[ Pn ]-v) , 



(29) 



E 



LSDA 



[PnA) 



jd 3 rp„(r)£ x L c SDA 



(W;r), (30) 



with £xc SDA the LSDA exchange and correlation energy 
density. The orbital densities entering in the SI term 
are those associated to the local orbitals <j>. As we have 
already mentioned, this functional needs to be minimized 
with respect to the </>'s, which are an implicit function of 
the total spin density p a . In the ASIC approximation 
these orbitals are not minimized, but taken as atomic 
functions. This means that in the present form the theory 
is not variational, in the sense that there is no functional 
related to the KS equation (|27fl by a variational principle. 
With this in mind we adopt the expression of equation 
(|28|l as a suitable energy, where the orbital densities are 
those given by the ASIC orbitals 



(31) 



In our implementation the LSDA KS energy E LSBA is 
directly available as calculated in the Siesta code [2^ and 
thus, only the second term of equation l(2*5|) needs to be 
calculated. This is easily done by calculating both U and 
£j C SDA on an atomic radial grid for each atomic orbital 
in the system. 
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D. ASIC and LDA+C/ 

We now compare our ASIC method with another 
atomic-like correction to LSDA, namely the LDA+C/ 
method HQ- In LDA+C/ one replaces the LSDA ex- 
change and correlation energy associated to the "corre- 
lated" orbitals (d or / shells), with the Hubbard-C/ en- 
ergy. Thus the functional becomes 

£ LDA+U [p(r)] = £ LSDA [p(r)] + £ U [{KJ] - £ DC [KJ] , 

(32) 

where the Hubbard energy _E U and the double counting 
term E BC depend on the orbital populations p^ of the 
correlated orbitals. Several forms for the LDA+C/ func- 
tional have been proposed to date. A particularly simple 
and transparent one 0, , which is also rotationally 
invariant, redefines the U parameter as an effective pa- 
rameter U c g = U — J and the functional takes the form 



E u -E 



DC 



^EE 

I ma 



t in in 



V v Ia P Ia 



(33) 



where in we separate out the index for the atomic position 
/ from the magnetic quantum number m, and introduce 
the off-diagonal populations p\° n = ~£ a /« (Cl^m*IO 
with P^* = Note that the LDA+C/ functional 

is SI free for those orbitals that are corrected. 

Although a rotationally invariant form of the ASIC 
potential can be easily derived, we assume here for sim- 
plicity that the system under consideration is rotationally 
invariant, or alternatively that we have carried out a ro- 
tation, which diagonalizes the p\^ n matrix. In this case 
the energy becomes simply 



E v -E 

with z4 CT = p\£ m 
potential 



DC 



2 



E & t 1 - & ] 



(34) 



It is then easy to compute the KS 



^LDA+U _ ^LSDA 



and the orbital energy 
r„ dE 



Im a 



77 - P 



I a 



P. 



/<I> 



dpi 



la LSDA , tj 



2"P B 



(35) 



(36) 



These need to be compared with the ASIC potential 
(equation l|l u pl and orbital energy 



Jo LSDA 



ap n 



;> , (37) 

where the last term follows from = G^T an< ^ f rom 
equation (|28l) . The main difference between the ASIC 
and LDA+C/ method is in the way in which unoccupied 
states are handled. In fact, while LDA+C/ corrects un- 
occupied states and pushes the orbital energies upwards 



by ~ C/ ff/2, ASIC operates only on occupied states, that 
are shifted towards lower energies by C%. This reflects 
the fact that the SIC is defined only for occupied KS or- 
bitals. An important consequence is that the opening of 
bandgaps in the electronic structures, one of the main 
features of both the LDA+C/ and ASIC schemes, is then 
driven by two different mechanisms. On the one hand in 
LDA+C/, gaps open up since occupied and unoccupied 
states are corrected in opposite directions leading to a 
gap of ~ Uqs- On the other hand ASIC is active only over 
occupied states and gaps open only if occupied and un- 
occupied bands have large differences in their projected 
atomic orbital content. Thus one should not expect any 
corrections for covalent materials where conduction and 
valence bands are simply bonding and antibonding states 
formed by the same atomic orbitals. This is for instance 
the case of Si and Ge. In contrast ASIC will be extremely 
effective for more ionic situations, where the orbital con- 
tents of conduction and valence bands are different. 

Finally, by comparing the corrections to the orbital 
energy of a fully occupied state, one finds 



U = 2a ($ 



V 



ml 



(38) 



which establishes an empirical relation between the Hub- 
bard energy and the ASIC correction. Since U is sensi- 
tive to the chemical environment due to screening [To| . 
while all the other quantities are uniquely defined by an 
atomic calculation, we can re-interpret the parameter a 
as empirically describing the screening from the chemical 
environment within the ASIC scheme. 



IV. RESULTS: EXTENDED SYSTEMS 

The test calculations that we present in this work are 
for two classes of materials: extended and finite. First we 
investigate how our implementation performs in the solid 
state. In particular we discuss the role of the parameter a 
in determining the bandstructure of several semiconduc- 
tors, considering both the KS band-gap and the position 
of bands associated with tightly bound electrons. 



A. Estimate of a for semiconductors 

The quasi-particle band gap E g in a semiconductor is 
defined as the difference between its ionization potential 
/ and electron affinity A. These can be rigorously calcu- 
lated from DFT as the HOMO energy respectively of the 
neutral and negatively charged systems. This actual gap 
cannot be directly compared with the KS band-gap E^ s , 
defined as the difference between the orbital energy of the 
HOMO and LUMO states of the neutral system. In fact, 
the presence of a derivative discontinuity in the DFT en- 
ergy as a function of the electron occupation establishes 
the following rigorous relation [20, 69] 



E n = E. 



KS 



(39) 
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with 



lim 



6E^ 



6n 



N+LU 



5n 



(40) 



N- 



This is valid even for the exact XC potential, and there- 
fore in principle one has to give up KS bandstructures as 
a tool for evaluating semiconductor band-gaps. The size 
of A xc is however not known for real extended systems 
and the question of whether most of the error in deter- 
mining Eg from Eg^ s is due to the approximation in the 
XC potential or due to the intrinsic A xc is a matter of 
debate. 

In general Si-free potentials bind more than LSDA and 
one expects larger gaps. Surprisingly, functionals based 
on exact exchange, prov ide KS gaps rather close to the 
experimental values |70l r71| . The reason for such a good 
agreement is not fully understood, but it is believed that 
the exact KS gaps should be smaller than the actual ones. 

With this in mind, we adopt a heuristic approach and 
we use the KS band-gap as a quality indicator for inter- 
preting the parameter a and for providing its numerical 
value for different classes of solids. Here we investigate 
the dependence of E^ over a and we determine the value 
for a yielding the experimental band-gap. Assuming that 
A xc does not vary considerably across the materials in- 
vestigated, this will allow us to relate a to the degree of 
localization in a semiconductor and to extract the value 
useful for ASIC to be an accurate single-particle theory. 

In figure n we present the band-gap of four representa- 
tive semiconductors as a function of a together with the 
value needed to reproduce the experimental band-gap. 
LSDA corresponds to a = and while a = 1 accounts for 
the full ASIC. In general Ey° increases as a increases, as 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

a 

FIG. 1: Calculated band-gap for NaCl, MgO,AlN and ZnO 
as a function of the parameter a. a — is the LSDA value 
and a = 1 accounts for complete atomic SIC. The lattice 
parameters used for the calculations are either the equilibrium 
LSDA or the experimental when available. 

a result of the stronger SIC. The Eg (a) curve is almost 
linear with a slope, which appears to be material-specific. 



For the most ionic compound, NaCl, the experimental 
gap is reproduced almost exactly by a = 1, i.e. by the 
full ASIC. This is somehow expected since the charge 
density of solid NaCl is rather close to a superposition 
of the Na + and Cl~ ionic charge densities. In this case 
of strongly localized charge densities the ASIC approx- 
imation is rather accurate yielding results substantially 
identical to those obtained with full self-consistent PZ- 
SIC |4l| . Indeed earlier calculations for LiCl demon- 
strate that the SIC band-structure is rather insensitive 
of the localized orbitals cf> once these have an atomic-like 
form. 

For the other compounds the localized orbitals </>'s are 
not necessarily atomic-like functions and deviations from 
a = 1 are expected. Interestingly we find that, for all 
the materials investigated, a value of around 1/2 repro- 
duces the experimental band-gap rather accurately. As 
an illustration, in table ^ we compare the experimental 
band-gap E^p to the calculated Ef s for ASIC (a = 1) 
and LDA, for several semiconductors ranging from ionic 
salts to wide-gap II- VI and III-V semiconductors. We 
also report the value of a = a* 



needed for E^ = Ef s . 



TABLE I: Experimental Ef p and KS Ef B band-gap (in eV) 



for a number of semiconductors 
both LSDA and ASIC (a = 1). 
port the value of a = a 



E^ s are calculated with 
In the last column we re- 



needed for E e * p 



E 



KS 



The 



lattice parameters used for the calculations are either the 
equilibrium LSDA or the experimental when available (in A). 
RS=rocksalt, WZ=wurtzite, ZB=zincblende. The value for 
the experiment al g aps are from the literature: a [7^|. b [73| . 

c m, d m, e Ea, / m, g m, h m, » m 



Solid 




Structure 


pcxp 


ttiKS-LSDA 
Eg 


ciKS-ASIC 


a 


LiCl 




RS (a =5.13) 


9.4 a 


6.23 


9.76 


0.89 


NaCl 




RS (a =5.63) 


8.6 6 


4.91 


8.51 


1.02 


KC1 




RS (a =6.24) 


8.5 C 


4.90 


8.51 


0.99 


MgO 




RS (a =4.19) 


7.8 d 


4.86 


9.36 


0.65 


CaO 




RS (a =4.74) 


7.08 d 


4.93 


9.28 


0.49 


SrO 




RS (a =5.03) 


5.89 e 


4.20 


7.80 


0.47 


A1N 


WZ 


(a=3.11, c=4.98) 


6.20 / 


4.47 


7.56 


0.56 


GaN 


WZ 


(a=3.16, c=5.13) 


3.39 9 


2.21 


5.03 


0.44 


InN 


WZ 


(o=3.54, c=5.70) 


0.7 h 


0.09 


2.09 


0.45 


ZnO 


WZ (a=3.23, c=5.19) 


3.43* 


0.85 


5.13 


0.57 


ZnS 




ZB (a =5.40) 


3.78* 


2.47 


4.90 


0.53 


ZnSe 




ZB (a =5.63) 


2.82* 


1.77 


3.53 


0.58 



Clearly for all the strongly ionic compounds (LiCl, 
NaCl and KC1) the full ASIC correction a = 1 repro- 
duces quite accurately the experimental gap and agrees 
with previous self-consistent SIC calculations 0|. For 
all the other compounds a value of around 1/2 is always 
adequate, confirming the initial choice of Filippetti and 
Spaldin. For these materilas we do not find any particular 
regularity In general a is large when the experimental 
gap is large, however there is no direct connection be- 
tween a and the ionicity or covalency of a compound. 
In fact, the improvement of the band-gap is not simply 
due to a rigid shift of the valence band, but usually cor- 
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responds to a general improvement of the whole quasi- 
particle spectrum. Examples for ZnO and GaN will be 
presented in the next section. 

As a further proof of this point in tabic ITT1 we present 
the valence band-width for the semiconductors investi- 
gated as calculated from LSDA A£ LSDA and ASIC for 
both a = 1 (A£ ASICl ) and a = a * (A£;ASic Q , y We also 

report the experimental values A_E° xp whenever avail- 
able, although a direct comparison with experiments is 
difficult, since these values are rather imprecise and some- 
times not known. The general feature is that ASIC pro- 
duces only minor corrections over LSDA, and that these 
corrections do not follow a generic trend. Thus, while 
for the nitrites ASIC always increases the band-width, it 
does just the opposite for KC1, SrO and CaO. 



TABLE II: Valence band experimental bandwidth A_EJ xp 
compared with those obtained from ASIC (a = 1) A_E ASICl , 
LSDA (AEv SDA ) and ASIC with the optimal a = a* from ta- 
ble [I] Ai? ASICQ * for a number of semiconductors. The lattice 
parameters used for the calculations are either the equilib- 
rium LSDA or the experimental when available (in A). The 
experimental values are from the literature (last column). 



Solid 


A££ xp 


A£ ASI Cl 


A£ LSDA 


A£ ASIC °* 


Reference 


LiCl 


4-5 


3.52 


3.06 


3.51 


[8JJ 


NaCl 


1.7-4.5 


2.11 


2.06 


2.11 


f81] 


KC1 


2.3-4.3 


1.09 


1.21 


1.09 


[81] 


MgO 


3.3-6.7 


5.16 


4.83 


5.06 


[82, 83] 


CaO 


0.9 


2.72 


2.89 


2.82 


[83] 


SrO 




2.21 


2.53 


2.39 




A1N 


6.0 


7.44 


6.27 


6.92 


[84] 


GaN 


7.4 


8.42 


7.33 


7.85 


[85] 
[86] 


InN 


6.0 


6.66 


6.01 


6.34 


ZnO 


~5 


5.66 


4.77 


5.54 


[87] 


ZnS 


5.5 


6.49 


5.57 


6.05 


[88] 
[88] 


ZnSe 


5.6 


7.14 


5.35 


6.38 


B. 


Wide- 


gap semiconductors: ZnO and GaN 



Having established a = 1/2 as an appropriate value 
for II- VI and III-V semiconductors, we now look at the 
whole band-structure (not just the fundamental gap) 
for a few test cases. Here we consider ZnO and GaN 
for which photo-emission data disagree quite remark- 
ably from LSDA calculations. In figure [2] we compare 
the band structure of wurtzite ZnO obtained respectively 
from LSDA and our ASIC. 

In ZnO, the valence band top (VBT) is predominantly 
oxygen 2p in character and the conduction band mini- 
mum (CBM) is essentially zinc 4s. With a value of ~ 0.5 
for the scaling parameter a, the ASIC band gap closely 
matches the experimental gap of E g =3.43 eV, whereas 
the LDA band-gap is very small (~ 0.85 eV). Some part 
of the LDA band gap error in ZnO can be traced to an 
underestimation of the semi-core Zn 3d states. The LDA 
binding energy for the Zn 3c? states is ~ 5.5eV while 
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FIG. 2: Calculated band structure of wurtzite ZnO obtained 
from LSDA and ASIC. Owing to the ionic character of ZnO 
each group of bands can be clearly labeled according to a 
single, dominant orbital character as shown. The VBT is 
aligned at eV. 



photoemission results place them at around ~ 7.8 eV. 
ASIC however rectifies the problem and is in very good 
agreement with experiment. This results furthermore in 
the removal of the spurious Zn3c;-02p band mixing seen 
in LDA. An additional feature is that the band-width of 
the valence band increases considerably as an effect of the 
downshift of the d manifold. Its worth mentioning that 
the positions of the Zn 3d levels obtained from ASIC in 
the case of ZnS, ZnSe, and ZnTe also agree remarkably 
well with experiment. 

The wide-gap III-V semiconductor GaN presents simi- 
lar phenomenology to that of ZnO. Figure[3]compares the 
band structure for wurtzite GaN obtained from LSDA 
and ASIC. When compared to X-ray photoemission spec- 
tra [83, the LSDA band structure of GaN has several 
shortcomings. Firstly, the band-gap between N 2p bands 
(VBT) and Ga 4s bands (CBM) is underestimated at 
around 2.2 eV against the experimental value of 3.4 eV. 
Secondly, the 3c? states of Ga are too shallow in LSDA, 
leading to a spurious 3d-2s hybridization. As a result 
the Ga 3d states overlap with and split the N 2s bands. 
ASIC rectifies the picture on both counts by improving 
the band gap and lowering the position of the Ga 3d 
bands with respect to the N 2s bands. 



C. Transition-metal oxide: MnO 

Transition metal oxides like MnO and NiO are char- 
acterized by partially filled 3d orbitals and an associated 
local magnetic structure. In particular the Mn 2+ ions 
in MnO are magnetic with a half-filled 3d shell. In the 
ground state, MnO is an A-type anti-ferromagnetic insu- 
lator in the intermediate charge-transfer Mott-Hubbard 
regime with a band-gap in the region of 3.8-4.2 eV. The 
VBT is expected to be of mixed Mn 3d-0 2p char- 
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FIG. 3: Calculated band structure of wurtzite GaN obtained 
from LDA and ASIC. The primary orbital character of the 
bands is indicated. The VBT is aligned at eV. 



actcr and the CBM pure Mn 3d in character. How- 
ever the LSDA description of MnO is flawed in sev- 
eral aspects most notably in describing MnO as a nar- 
row gap (E g = 0.92 eV) Mott-Hubbard insulator with 
both the VBT and CBM composed of purely of Mn 3d 
states. This is due to the severe underestimation of d 
electron binding-energies in LSDA. The calculated anti- 
ferromagnetic band-structures of MnO from LSDA and 
ASIC (a — 1/2) are shown in figure 0] 




FIG. 4: Calculated band structure of anti-ferromagnetic MnO 
obtained from LSDA and ASIC. In our calculation we obtain 
an LSDA bandgap of ~0.65 eV whereas the ASIC bandgap is 
much improved at ~3.5 eV. The VBT is aligned at eV. 

Note that these are for the rhombohedral unitcell 
with 4 atoms per cell. The two Mn ions are anti- 
ferromagnetically aligned and the oxygen ions are non- 
magnetic. This results in a layered ferromagnetic order 
of the (111) planes, which in turn are anti-ferromagnetic 
coupled to each other. Also in this case, ASIC is a con- 
siderable improvement over LSDA. The size of the fun- 



damental gap now resembles the experimental one and 
the VBT recovers some p character. 



RESULTS: MOLECULES 



A. Ionization potentials 

In view of the fact that the ASIC method gives im- 
proved eigenvalue spectra for several solid state systems, 
it is worth taking a cautious look at how it performs with 
molecules. This is particularly important for assessing 
whether the ASIC scheme can be adapted to work in 
DFT electron transport schemes based on the KS spec- 
tra [HiiJ. In exact KS DFT only the highest occupied 
orbital eigenvalue ( e HOMO ) has a rigorous physical in- 
terpretation and corresponds to the negative of the first 
ionization potential [l3 . |2C| . More generally, for a N 
electron system, the following equations hold in exact 
KS-DFT ' 



HOMO 



HOMO 



(M) = -7 N for (N — 1 < M < N) (41) 



(M) = -A N for {N < M < N + 1) (42) 

where —In and —An are the ionization potential (IP) and 
the electron affinity (EA) respectively. Therefore we start 
our analysis by looking at these quantities as calculated 
by ASIC. Also in this case we investigate different values 
of a. However here we limit ourselves only to a = 1 
(ASICi) and a = 1/2 (ASIC 1/2 ). 

In table ITTT1 and figure [3] we compare the experimental 
negative IP for several molecules with the corresponding 
(e HOMO ) obtained using LSDA and ASIC. It is clear that 
LSDA largely underestimates the removal energies in all 
the cases and that the values obtained from ASIC1/2 
are also consistently lower than the experimental value. 
However, as made evident by the figure the agreement 
between ASICi and experiments is surprisingly good. In 
fact the mean deviation 8(X) {X = LSDA, ASIC1/2, 
ASICi) from experiment 



8{X) = \ 



E 



N 



HOMO.i 

-X 



Expt 



N 



is 3.56 eV for LSDA, 1.69 eV for ASIC 1/2 and only 0.58 
eV for ASICi (N runs over the molecules of table ITTTfl . 
It is worth noting that we have excellent agreement over 
the whole range of molecules investigated going from N2 
to large fullerenes Ceo and C70. 

For comparison in figure [5] we have also included re- 
sults obtained with a full self-consistent PZ-SIC approach 
|35| . Surprisingly our atomic approximation seems to 
produce a better agreement with experiments than the 
self-consistent scheme, which generally overcorrects the 
energy levels. This is a rather general feature of the PZ- 
SIC scheme and it is generally acknowledged that some 
re-scaling procedure is needed [9^, . 
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TABLE III: Experimental Ionization potential (IP) com- 
pared to calculated HOMO eigenvalues for neutral molecules. 
Columns 3 and 4 present the results from ASIC with respec- 
tively a = 1/2 and a = 1. The experimental data are taken 
from reference |91|. 



Molecule 




E HOMO (eV) 




-IPfeV) 




LSDA 


ASICl /o 


ASICi 


Exp eriment 


CH 3 


-4.65 


-7.34 


-10.06 


-9.84 


NH 3 


-5.74 


-8.21 


-10.79 


-10.07 


SiH 4 


-7.95 


-10.14 


-12.41 


-11.00 




-6.28 


-8.00 


-9.74 


-10.51 


SiCH 4 


-5.89 


-7.57 


-9.35 


-9.00 


CH3CHCI2 


-7.23 


-8.97 


-10.72 


-11.04 


C4H4S 


-5.95 


-7.65 


-9.35 


-8.87 


C2H6S2 


-5.56 


-7.54 


-9.53 


-9.30 


Pyridine 


-4.83 


-6.57 


-8.31 


-9.60 


Benzene 


-5.92 


-7.59 


-9.28 


-9.24 


Iso-butene 


-5.39 


-6.98 


-8.6 


-9.22 


Nitrobenzene 


-6.49 


-8.76 


-10.67 


-9.92 


Naphthalene 


-5.49 


-7.04 


-8.59 


-8.14 


Ceo 


-5.06 


-6.53 


-8.02 


-7.57 


C70 


-4.92 


-6.40 


-7.89 


-7.36 




FIG. 5: Experimental negative ionization potential IP com- 
pared to the calculated HOMO eigenvalues for molecules. The 
experimental data are from reference |9l|. while the star sym- 
bol represents full PZ-SIC calculations from reference |3a| . 



B. Electron affinities 

In Hartree Fock theory where Koopmans' theorem 
holds the lowest unoccupied molecular orbital 

(LUMO) energy (e LUMO ), corresponds to the vertical EA 
of the N electron system, if one neglects electronic relax- 
ation. No such interpretation exists for ( e LUMO ) m DFT 
and so the EA is not directly accessible from the ground 
state spectrum of the N electron system. However, as 
equation l|42|) indicates, the EA is in principle accessi- 



ble from the ground state spectrum of the N + 1 — / 
(0 < / < 1) electron system and asserts in particular 
that it must be relaxation free through non-integer oc- 
cupation. Unfortunately, the LSDA/GGA approximate 
functionals usually perform rather poorly in this regard 
as the N + 1 electron state is unbound with a positive 
eigenvalue. So one resorts instead to extracting electron 
affinities from more accurate total energy differences |95j , 
or by extrapolatin g th em from LSDA calculations for the 
N electron system |9(| • This failing of approximate func- 
tionals has been traced in most part to the SI error and 
so SIC schemes are expected to be more successful in 
describing the N + 1 electron spectrum. 

In table ITVl we compare HOMO energies (denoted as 
of several singly negatively charged molecules 
with the experimental electron affinities. We also re- 
port the LUMO energies for the corresponding neutral 
species (denoted as e^ UMO ). LSDA relaxed geometries 
for the neutral molecule are used for both the neutral 
and charged cases. We find that various En+i 10 obtained 
from ASICi once again are in reasonably good agreement 
with corresponding experimental electron affinities while 
LSDA and ASIC1/2 continue to be poor even in this re- 
gard. In this case 5(X) stands at 4.1 eV , 2.31 eV and 
0.83 eV for LSDA, ASIC1/2 and ASICi respectively. No- 
tice that e^MO from LSDA is positive in most cases as 
the states are unbound. 

In figureElwe present our data together with £^+5^° as 
calculated using the PZ-SIC j3^. Again ASICi performs 
better than PZ-SIC, that also for the EA systematically 
overcorrects. 




FIG. 6: Experimental negative electron affinities (-EA) com- 
pared to calculated HOMO eigenvalues of negative radicals. 



C. Vertical excitations 

Having shown that ASIC offers a good description of 
both IP and EA for a broad range of molecules, we turn 
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TABLE IV: Calculated HOMO eigenvalues for singly negatively charged molecules compared to experimental negative electron 
affinities (-EA). Columns 6,7 and 8 present the LUMO eigenvalues for the corresponding neutral species. 

Molecule eggf°(cV) Exp. -EA (eV) e^ UMO (eV) 





LSDA 


ASIC 1/2 


ASICi 




LSDA 


ASIC 1/2 


ASICi 


CN" 


0.84 


-0.79 


-2.48 


-3.86 


-8.13 


-9.03 


-9.42 


C 2 H" 


0.94 


-0.80 


-2.68 


-2.97 


-6.91 


-7.38 


-7.48 


CH 3 S" 


2.42 


0.65 


-1.14 


-1.87 


-5.20 


-5.31 


-5.34 


Uri 


000 
3.s2 


1 no 

i.uy 


1 on 
-1.80 


-1.83 


Old 

-U.16 


-0.4,5 


-U.69 


S1H3 


4.61 


3.13 


1.61 


-1.41 


-2.66 


-3.30 


-4.07 


HOO" 


3.10 


-0.07 


-3.34 


-1.08 


-5.30 


-6.14 


-6.40 


NH 2 " 


3.83 


1.51 


-0.98 


-0.77 


-5.27 


-4.80 


-4.39 


CH2 


3.07 


1.21 


-0.45 


-0.65 


-3.80 


-3.84 


-3.91 


CH3CO- 


2.90 


1.76 


0.40 


-0.42 


-2.94 


-3.88 


-4.85 


CHO" 


3.55 


2.02 


0.42 


-0.31 


-3.30 


-4.40 


-5.51 


CH3 


4.15 


1.99 


-0.34 


-0.08 


-2.73 


-2.59 


-2.47 


^60 


0.03 


-1.19 


-2.45 


-2.65 


-3.44 


-4.66 


-5.90 


C70 


0.00 


-1.22 


-2.47 


-2.73 


-3.17 


-4.41 


-5.66 



TABLE V: Orbital energies of N2 calculated with various 
methods. The results for Hartree-Fock and SIC are from ref- 
erence |3lT [. Experimental results are from reference |9?| . 



Orbital 


HF 


SIC 


D-SIC ASICi 


ASICi /a 


LSDA 


Exp. 


2a g 


-41.49 


-38.86 


-37.85 


-38.29 


-33.22 


-28.16 




2a u 


-21.09 


-20.27 


-16.44 


-18.42 


-15.64 


-12.93 


-18.75 


3(Jg 


-17.17 


-17.39 


-13.88 


-14.01 


-11.70 


-9.90 


-15.58 




-16.98 


-16.33 


-16.68 


-15.97 


-13.74 


-11.54 


-16.93 



our attention to the remaining vertical ionization poten- 
tials. As mentioned before, KS-DFT lacks of Koopmans 
theorem, and therefore the KS energies are not expected 
to be close to the negative of the vertical ionization po- 
tentials. However, at least for atoms, the introduction 
of SIC brings a remarkable cancellation between the neg- 
ative relaxation energy and the positive non-Koopmans 
corrections [l^. For this reason the SIC KS eigenval- 
ues are a good approximation to the relaxed excitation 
energies. As an example, in table E] we present the or- 
bital energies calculated with ASICi and ASIC1/2 for 
the No molecule. These are compared with experimental 
data [93 and orbital energies obtained respectively with 
Hartree-Fock (HF), self-consistent SIC, and SIC where 
molecular orbitals are used instead of localized orbitals 
(D-SIC) H3. 

Remarkably ASICi seems to offer good agreement over 
the whole spectrum, improving considerably over LSDA 
and in some cases even over SIC and HF results. This im- 
provement is not just quantitative, but also qualitative. 
For instance while rectifying the LSDA spectrum of the 
N2 molecule, ASICi preserves the correct order between 
3cr g and lir u orbitals, which are erroneously inverted by 
both SIC and HF. So why does ASIC perform better than 
the other methods with regards to removal energies? In 
LSDA, electron relaxation typically cancels only half of 
the non-Koopmans contributions, resulting in energies 



TABLE VI: Orbital energies for CO, HF and H 2 calculated 
with LSDA and ASICi. The experimental results are from 
reference |98| and references therein. 



Molecule 


Orbital 


LSDA 


ASICi 


Exp. 


CO 


5a 


-8.74 


-12.85 


-14.01 




In 


-11.54 


-16.64 


-16.91 




ia 


-13.97 


-19.36 


-19.72 


HF 


In 


-9.83 


-16.96 


-16.19 




•V 


-13.61 


-19.68 


-19.90 


H 2 


lbi 


-7.32 


-13.38 


-12.62 




3ai 


-9.32 


-14.66 


-14.74 




lb 2 


-13.33 


-18.03 


-18.55 



that are too shallow [16j. In contrast HF lacks energy 
relaxation and the orbital energies are too deep. The 
reason why ASICi performs better than self-consistent 
SIC is less clear. As a general consideration, also for 
the case of vertical ionization energies self-consistent SIC 
seems to overcorrect the actual values. Thus the SIC 
potential appears too deep, and the averaging procedure 
behind the ASIC approximation is likely to make it more 
shallow. 

As a further test we calculated the orbital energies for a 
few other molecules and compared them both with LSDA 
and experiment [98| . These are presented in table IVII 
Again the ASICi results compare rather well with ex- 
periment, and we can conclude that the ASIC method 
offers a rather efficient and inexpensive theory for single 
particle vertical excitations. 



D. HOMO-LUMO gap and discontinuity of the 
exchange and correlation potential 

We are now in a position to discuss the HOMO-LUMO 
gap in ASIC. As already mentioned, even for the ex- 
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^LUMO 



r HOMO 



act XC functional, the KS gap E^ s = e"^"^ — e 
does not account for the actual quasi-particle gap E g = 
In — j4n- This in turn is the sum of Eg and the dis- 
continuity of the exchange and correlation potential A xc . 
Equivalently 



A xr = lim gHOMO ,LUMO 



-N+f 



■ N 



(43) 



i.e. A xc is the discontinuity in the eigenvalue of the 
LUMO state at N. Therefore, in order to extract the 
actual gap from the KS gap, provided that the spectrum 
is reasonably well described at integer electron numbers 
N, what remains is to model the derivative discontinuity 
at TV and ensure that e^.y?° is relaxation free for (0 < 
/ < 1). Local and semi-local (LSDA/GGA) XC func- 
tional lack such a discontinuity, while self-interaction 
corrections are able to restore it, at least in part. For 
instance the PZ-SIC scheme is successful in this regard. 




0.4 0.6 

p(\e\) 

FIG. 7: Ionization curve for the ethylene (C2H4) molecule as 
the occupation of the HOMO state is varied from to 1 in 
going from the ionized C2HJ to the netural C2H4 state. 

In figure \7\ we illustrate the ionization curve for 
the ethylene (C2H4) molecule as the occupation of the 
HOMO state is varied from to 1 in going from the ion- 
ized C2H4 to the netural C2H4 configuration. It is seen 
that among the three schemes presented, only the PZ-SIC 
scheme approximately models the behaviour required by 
the equation igTJ. The ASIC HOMO eigenvalue roughly 
agrees with the PZ-SIC eigenvalue at integer occupation 
but behaves linearly through non-integer values. Thus we 
find that the derivative discontinuity for the molecule is 
smoothed out in ASIC, which still connects continuously 
different integer occupations. This is one of the limita- 
tions of the atomic representation employed in ASIC. 

In view of the foregoing discussion, the actual size 
of the HOMO-LUMO gap in ASIC becomes significant 
with a direct bearing on the physics described. Ide- 
ally, we want e LUMO (jV) (LUMO for the iV-electron sys- 
tem) to be as close to e HOMO (N + 1) so that the range 
of eigenvalue relaxation through fractional occupation 



TABLE VII: HOMO-LUMO gap obtained from ASIC com- 
pared to the LSDA value. The values marked with * corre- 
spond to unbound LUMO levels. 



Molecule 



(eV) 





LSDA 


AST(~!, ,., 

-ri-Olv^i/2 


ASICi 


v 1 1 ; 


1.92 


A 7^ 


7 to 


NHq 


7.1* 


9.29* 


11.61* 


S1H4 


8.44* 


9.68 


10.94 


C2H4 


5.81 


6.59 


7.38 


SiCH 4 


6.19 


7.07 


8.06 


CH3CHC12 


5.79 


6.84 


7.88 


C4H4S 


4.46 


5.13 


5.8 


C2H6S2 


4.44 


6.02 


7.6 


Pyridine 


3.85 


4.56 


5.26 


Benzene 


5.22 


5.9 


6.59 


Iso-butene 


4.88 


5.56 


6.26 


Nitrobenzene 


3.25 


4.03 


4.42 


Naphthalene 


3.36 


3.83 


4.29 


Ceo 


1.62 


1.87 


2.12 


C70 


1.75 


1.99 


2.23 



numbers M G {N,N + 1) is minimized. Looking at 
columns 6,7 and 8 in table ITVl however, we see that for 
almost all the molecules, this is hardly the case. The 
agreement between e LUMO (]\r) and -EA from experiment 



c HOMO 



(N + 1)) for ASICi is quite poor implying 



a considerable energy range spanning fractional particle 
number. We still expect this energy range to be smaller 
for ASICi than LSDA. It is also apparent from the table 
El that (AO usually differs from e£g^°(iV) and 

in fact by considerable magnitudes in some cases. Thus 
the ASICi "correction" to the empty LUMO state does 
not vanish in contrast to the PZ-SIC scheme where, by 
definition, the empty eigenstates are SIC free. 

Since the SIC operator ^Xsic ^ s constructed in an 
atomic orbital representation, the correction to any KS 
eigenstate -0n either filled or empty 



tE$&c = «KsiclO 



(44) 



is not necessarily zero unless tp^ only projects onto empty 
atomic orbitals. Also this correction to the LUMO with 
respect to the LSDA is negative in most cases, exceptions 
being NH 2 and CH 3 where it is desirably positive. Thus 
the fundamental HOMO-LUMO gap in ASIC is a combi- 
nation of both the HOMO and LUMO corrections. Table 
IV 111 shows how this combination works out in ASIC1/2 
and ASICi when compared to LSDA. The molecular test 
set is the same as that in table Hill 

We see in almost all cases the ASIC gap is system- 
atically larger than the LSDA one. This is expected 
because the correction to the HOMO is usually much 
stronger than that to the LUMO. In general, ASIC is 
expected to work well for systems where the occupied 
and un-occupied KS eigenstates of the extended system 
have markedly different atomic orbital signatures being 
derived predominantly from filled and empty atomic 



14 



orbitals respectively. In such a case, the ASIC correction 
to the empty states would be nullified in being scaled 
by near-zero atomic orbital populations. In some cases, 
provided phase factors combine suitably, the correction 
to the empty states can even be positive with respect to 
the same in LSDA. 



E. Final Remarks 

Before we conclude, we discuss some general properties 
of the ASIC method which are relevant to any orbital de- 
pendent SIC implementation and also some possible pit- 
falls. As with other SIC schemes, ASIC is not invariant 
under unitary transformations of the orbitals used in con- 
structing the SIC potential. Thus the ASIC correction is 
likely to change as the atomic orbitals used for project- 
ing onto the KS eigenstates of the system are rotated 
or transformed otherwise. Unlike the Perdew-Zunger 
method however, there can be no variational principle 
over all possible unitary transformations of the atomic 
orbitals because in the general case they do not repre- 
sent the Hamiltonian of the system under consideration. 
This also implies that if the scheme is used with a system 
that is already well described by LSDA, the "correction" 
additional to the LSDA result does not necessarily van- 
ish. Simple metals and narrow gap systems are likely 
candidates for this scenario. 

Furthermore, its pertinent to mention that ASIC be- 
comes ineffective if not counterproductive for materials 
with homonuclear bonding, in which valence and conduc- 
tion bands have the same atomic orbital character. In 
this situation the ASIC potential will shift the bands in 
an almost identical way, without producing any quantita- 
tive changes, such as the opening up of the KS gap. Note 
that this is a pitfall of the ASIC approximation, which 
distinguishes occupied from empty states only through 
their projected atomic orbital occupation, but not of the 
SIC in general. Typical cases are those of Si and Ge. 
The KS gap in Si goes from 0.48 cV in LSDA to only 
0.57 eV for ASIC1/2, while Ge is a metal in both cases. 
In addition the LSDA calculated valence bandwidths of 
12.2 cV for Si and 12.8 eV for Ge, in good agreement 
with experiments, are erroneously broadened to 14.3 eV 
and 14.8 eV respectively. 



VI. CONCLUSIONS 

In conclusion, we have implemented the ASIC scheme 
proposed by Filipetti and Spaldin within the pseudopo- 
tential and localized orbital framework of the Siesta code. 
We have then investigated a broad range of semiconduc- 
tors and molecules, with the aim of providing a reason- 
able estimate for the scaling parameter a. We found that 
a = 1, which accounts for the full atomic SI, describes 
surprisingly well ionic semiconductors and molecules. In 
particular for molecules, both the IP and the EA can be 
obtained with good accuracy from the HOMO KS eigen- 
values respectively for the neutral and singly charged 
molecule. This makes the ASIC scheme particularly 
suited for application such as quantum transport, where 
the position of the HOMO level determines most of the 
I-V curve. 

In contrast III-V and II- VI semiconductors are bet- 
ter described by a — 1/2, which corrects the atomic SI 
for screening. This makes ASIC1/2 an interesting effec- 
tive band theory for semiconductors. The relation of the 
present scheme with the fully self-consistent SIC methods 
has been emphasized, and so has been that with LDA+J7. 
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